NOTICE 


THIS DOCUMENT HAS BEEN REPRODUCED FROM 
MICROFICHE. ALTHOUGH IT IS RECOGNIZED THAT 
CERTAIN PORTIONS ARE ILLEGIBLE, IT IS BEING RELEASED 
IN THE INTEREST OF MAKING AVAILABLE AS MUCH 
INFORMATION AS POSSIBLE 









Dec. 5, 1980 


Planetary Science Institute 
Science Applications, Inc. 


283 S. Lake Ave. 
Pasadena , / CA 



i 


% 
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This is the Pinal Report on Tasks 3 and 4 of tJie program 
"Photogeological Constraints on Lunar and l^lanetary Volcanism” 
which was performed between Dec. 7, 1978 and Dec. 7, 1979. 

The objectives of the program were to develop' an understanding 

of the physical mechanisp of thermal erosion by geophysical 

w" 

fluids. 

In Task 3A we attempted to develop a Conceptual under- 
standing of thermal erosion. We developed a simple one dimen- 
sional model of the process for a substrate with temperature 
^dependent viscosity in which erosion was specified to occur 
'when the substrate reached a critical value. This model pre- 
dicted a square root law dependence of erosion on time whereas 
common sense considerations indicate a roughly linear relation 
ship. The model was modified to include a convection term 
yielding the expected linear relationship between erosion and 
time, 

I 

In Task 3 b we explored some simple two dimensional models 
in which the downstream dimension was explicitly included in 
the model. All the models examined had simple analytical 
solutions and were primarily developed to improve our under- 
stranding of . different aspects of the physics of the thermal 

j; 

erosion process, in the first 2D model we calculated the 
thermal and velocity fields for a steady state flow of a fluid 
of temperature - independent viscosity with a constant heat 



influx at the base of the flow^ The calculated veldclty has 
a parabolic dependence on depth. However, the temperature 
depends linearly on downstream distance with both quartic 
and quadratic components in depth dependence. We also ex- 
amined thd Conditions for turbulent and fully developed flow 
in the context of this model. It appears that lava erosion 
flows may be turbulent of leuninar but will almost certainly 
be fully developed. 

Our second simple two dimensional model included an 
examination of the effect of a finite yield stress on flow. 
Flowing lava behaves as if it haS a finite yield stress forming 
'labs of finite thickness. We showed schematically that thermal 
erosion by a Bingham fluid results in initial accretion (de- 
position) of material; net erosion will only occur with a 
sustained flow. We regard this concept as important to 
understanding the reasons for lava accretion and erosion in 
geologic settings. 

in a final 2D model we attempted to account in a pure.’y 
phenomenological way for the effects of turbulence and mech- 
anical erosion in a flow. We adopted criteria for mechanical 
erosion related to the basal stress and we parameterized the 
effect of entrained eroded material on fluid properties. A 
predibted downstream slope profile W5is generated for steady 
state erosion. 
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In Task 4 ye extended the earlieir work to include labora- 
tory simulation; and two dimensional finite element models. 

In the laboratory simulations with hot wax (Task 4A) we con- 
firmed that when a hot fluid crosses an erodible substrate, 
accretion occurs initially and erosion will only occur if the 
flow is sustained. We also showed that the amount of erosion 
is greates% nearest the source and the crossover point between 
net erosion and accretion propagates downstream with time. 

The evolution of the cross section of a thermal erosion channel 
was also examined. This work was performed by Jose Helu, 

In Task 4B we developed two dimensional finite element models 
to simulate lava erision. This work was initially carried out 
by Wm. James Roberts with the assistance of Jose Helu. 

However, when Roberts left the Institute in September of 1979, 
the work was continued under the direction of the Principal 
Investigator with the assistance and support of Stephen J. 
Keihm. 

Model results were generated for lava flows with three 
different velocities (lOcm/sec, lOOcm/sec, and lOJOcm/sec) 
each lasting for 96.7 hours. Plots were generated at inter- 
vals of approximately 20 hours for each model and showed the 
propagation of isotherms relative to the original interface 
between a lava flow and substrate. Plots of single models 
showed how the position of the 600°c isotherm varied with time 
in a given model and contrasted the loci of the 60 0°C isotherm 
at the end of each of the three experiments. 
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Xn general as the velocity o£ the flow was increased 
the isotherms propagated more rapidly into the substrate which 
equates with more rapid erosion. Propagation is most rapid 
near the source of fluid. A logical next phase in these 
modelling efforts is to generate an explicit description for 
removal of material which could be ^bmpared with the experimental 
results but the best method for doi^fig tills is by no means 
clear. 

Data and concepts presented in this report are being in- 
cluded in an article on thermal erosion of martian channels 
which is in preparation and will be submitted to Icarus . 


]1 




1 


TABLE OP CONTENTS 


PART I 

;7 

Central Volcanic Constructs as Constraints on 
the Thermal Evolution and Regional Tectonic^ 
of the Moon and Terrestrial Planets 

Origins of S\ib-kli6meter Lunar Craters: Im- 

plications for Mare Basalt Petrogenesls 

APPENDICES TO TASKS 1 & 2 

PART II 

TASK 3: Lava Erosion: Theoretical Concepts and Simple 

Analytical Models 

TASK 4 ; Lava Erosion : Laboratory Simulation and Two 

Dimensional Finite Element Models 

APPENDICES TO TASKS 3 & 4 


TASK 1: 

1 . 

TASK 2: 






TASK 3: THEORETICAL MODELLING STUDY OF LAVA EROSION 

The use of several simple theoretical modelling approaches 
in simulating lava erosion were evaluated under this task. 

At the outset, when we planned this program, we perceived 
that a one dimensional model, with properties invariant in 
the downstream direction, was not suited to modelling a pro- 
cess with inherent variations in ^phe direction of flow. 

However, such an approach had been taken by previous workers 
and we decided to try it in order to establish a link with 
previous work. The one dimensional studies turned out to be 
limited in scope; we have also examined some highly simpli- 
fied two dimensional models as a part of this task. The re- 
sults of these investigations are summarized below. 

A. simple One Dimensional Model 

In the treatment of lava erosion given by Carr 1974, the 
complex intertwined relationships. of the thermal and velocity 
fields in the downstream dimensions are ignored. The lava 
erosion pfoblem reduces to a simple problem of thermal dif- 
fusion. The problem was examined for a channel with an 
initially rectangular cross-section. 

Carr originally treated the lava erosion process in terms 
of a yield stress. When the sxibstrate is heated to the yield 
teraperature the material of the channel walls immediately 
participates in the flow. The thermal burden placed on the 
flow by eroded material was not assessed; neither was the 
rate at which material could be assimilated and mixed. 
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Both of these effects are very difficult to model using the 
yield stress formulation. For this reason, we first looked 
at a one dimensional erosion model in which substrate changes 
from solid to liquid over a range of temperatures and not 
at a single yield point. 

1) simple one dimensional numerical modelling ■■ temperature 
dependent viscosity - conduction only 
This model was originally developed as part of a pre- 
posal pilot study. The model is described in detail in 
Appendix 3.1 and the results illustrated in Fig. 3.1. 

The depth of erosion, was (somewhat; arbitrarily) defined 
as the depth at which the substrate velocity falls below 
some critical value. This curious definition is needed be- 
cause the semi-infinite character of the boundary conditions 

permits no net removal or addition of substrate material. 

The model probably overestimates the time needed to 

remove one meter of substrate as it incorporates no mechanism 
for bringing fresh hot lava into contact with the substrate. 

The amount of erosion is proportional to the square root 
of time over a broad range of times and lava temperatures. 

This behavior is to be expected from the constraints of the 
model as described above but is not expected for real world 
lava erosion which is likely to more closely approach a linear 
law. Downstream variation in erosion and accretion is not 
part of the model. Finally, we did not include the temperature 
dependence of rock conductivity (Appendix 3,2), 








TIME (SECONDS) 


Pig# 3.1 Results of simp3;!& temperature dependent viscosity 
model of lava erosion. 










?) simple one-»dlmensional nuroerical modelling - temperature- 
dependent viscosity - conduction and convection 
Some of the deficiencies of the "one-dimensional-conduction- 
only" model can be rectified by including a cpRvectiop term to 
permit faster transfer of heat than allowed by conduction alone. 

We have assumed that material at all levels above the level with 
temp. and therefore,; viscosiuy are mixed in a time 

Interval t^. ''Our specific implementation of this model also allows 
the addition of hot lava from upstream. 

A grid of 14 nodes (depths) at which temperature is defined 

(i 

\ 

is set up. Another similar grid is set up "upstream" of it and 
filled with hot fluid. At each time step, first the convection 
and conduction terms are evaluated and new temperatures are de- 
fined, all using only the main gs?id. Next, hot fluid is intro- 
duced from the upstream "reservoir” grid, as a function of the 
temperature in each main grid element. This simulates in a crude 
fashion the viscosity-dependent addition of hot lava from upstream. 
If a node is below no addition of hot lava takes place. 

If a node is above then the node is completely replaced 

by hot lava. For '^thr^^*^*^hiv' amount of replacement of 

partially melted material with hot lava depends linearly on T 
(Fig. 3-2) . New main grid temperatures are then defined taking 
this into account, and the cycle repeats. 

In Fig. 3.3 and 3.4 we show some results calculated with this 


model. Because the parameters of hot wax were used in this calculation 
the results cannot be directly compared with results shown in 
Fig. 3.1. However, the total erosion depends linearly on time ^ 


Temperature 


Fig. 3 


,2 The percent replacement by hot lava as a functio^^^^^ 
of lava temperature for the one dimensional mode 
with conduction and convection. 
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Fig. 3.3 Erogion as a Function of Time for One Dimensional 
Model with Gortduction. The Result is Linear as 
Expected . 









Fig. 3.4 Temperature Profile as a Function of Time for Simple 
One Dimensional Model with Conduction and Convection 




as expected in contrast ylth the simple conduction model 
(Section 3hl) i fhe typical temperature profile using the 
hot wax parauneters (Fig. 3.4) shWs a sharp break in the 
temperature gradient at the fluid base. 

The methods used above to simulate non-conductive heat 
transport are obviously crude. The mixing model could be 
improved upon using numerical or parameterized theories of con- 
vection but the simulation of advection of heat from upstreeum 
is obviously not capable of rigorous treatment in purely one 
dimensional terms. The two dimensional treatment is needed 
to accomplish this. 

B. Simple Two Dimensional Models 

1) Simple two dimensional models - laminar flow 

\ \ . 

Simple models with analytic solutions for the thermal 

(I - 

and velocity fields are sometimes useful for understanding 
the physics of a flow process. Such a model has these at- 
tributes: the flow is assumed to be in a steady state and 
basal heat flow is uniform at all places along the flow. 

We also ignore the heat generated by fluid friction and 
assume that fluid viscosity is independent of temperature. 

The geometry of this model is shown in Pig. 3.5(a). 

Velocity Field 

The Navier Stokes equations describing the flow of mass 
and momentum in an incompressible fluid with constant viscosity 
are given in Bay iey Owen and' Turner (Eq. 5.15) * 


SJtfc.Jfe'tJUj.ii.'' 



Fig. 3.5a Geometry o£ simple 2 two dimensional laminar 
flow model o£ fluid viscosity indep of 
temperature. 
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where p « density \ 

x,y a positional coordinates 
u,v * velocities 


■ component of body force along 


p » pressure 
p 3 dynamic viscosity 

Under conditions of steady flow along a plane: 

0 V * 0 everywhere 

u “ u(y) and is therefore independent of x 


3u 


=0 everywhere 


~ = is independent of x 
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3.B.1 


3.B.2 


3.B.3 

3.B.4 

3.B.5 


3.B.6 

3.B.7 


= is independent of x 
oy 

F = pg sin(j>; F^ » pg cos, 

X y 

Substituting these conditions in C3.B.1) ' and C3.B,2) we obtain 


3.B.8 

3.B.9 


3^u . X 

y_ — r » pg sin(p 

3x 








3y 


pg cos 
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3.B.10 

3.B.11 
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Integrating equation (3.B.10) once 


y + c. 


3.B.12 


At y • 0 1 ^ • 0 f therefore C * 0 


Integrating again 


(yn^ - y ) 


3.B.14 


The mean velocity 


/ “ u 


dy 


pg sing 




pg sind> yr 


» j 0 max 


3.B.15 


Thermal Field 


For an incompressible fluid the energy equation for the 


flow can be expressed using equation 5.42 (b) of Bay ley Owen 


and Turner as 


i - " ’ w ' ”"^7 ^ 


3.B.16 


The energy equation shows that convective heat transfer due 


to fluid motion is balanced by work due to volumetric changes, 


the heat conducted through the fluid and viscous dissipation. 


It doesn't contain potential energy terms as these are implicit 


within the other energy terms . 
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Under the conditions of steady flow down a plain then 
V “ 0 everywhere. Other conditions are that the basal heat 

flow (q_) is constant and <|) is negligible. We adopt a 

“ 2 

a 

further contraint on a solution that j « 0 throughout 
the flow which is the condition for a steady state fully- 


developed flow. Equation (3.$:. 16) reduces to: 
3T 


u 


3X 


3^ 
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3. 


Where aa 


R 


pc. 


a thermal diffusivity with boudary conditions: 




‘ 3 y' 


(||)y = 0 a 0 


From the velocity solution equation (3.B.14) 


U - 


2 . 2v _ „ 2 „2, 


0 ” y V ' ''1 ‘^0 

We will look for a solution satisfying (3.B.14, 3.B.17, 

3. B. 18), and 3.B.19 
T a T^G(y) + T2X 

where <^(y = yg) * 0, such that the surface temperature of 
the lava equals the temperature of the plane on which it 
flows. 


(y« - y") 


II = T • 

3X ^2' 


a T, 


3^T 


Evaluating derivatives of (3. B. 20) 
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3y 

Substituting in equation (3;B.17) 
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ay - -1= ’’1 


P^^(y) 


UT2 = 


(y) 
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G (y) = 


aT, 


(yn - y ) 
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B.17 

B.18 

B.19 

B.14A 

B.20 


B.21 

B.22 
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Integrating equation (3.B.22) once: 
sNy) ■ Cl (y,2y - yVj) * 

At y » 0 {■—) * 0 and, substituting in (16) ■ 0 

3T 
3y 

ct3 ^s 
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At y ■ Yq (|^) » "*33/^ substituting in (3.B.23) 


2R yQ- 


3.B. 


3.B. 


Now integrating (3, B. 23) with set to zero. 

s(y) - si • Ci(y,^ h ~h* °2^ 

and substituting for T 2 


G(y) 


3ag 


2H y^ aT^ 


- fe * °2> 


at y » 0, G(y) = 0 and D 2 * “^ 2 ' - substituting in 

equation (13) for T^, T 2 and G(y) 


3.B 


3.B 


T 
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pg sin(()J 
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The computer program in Appendix 3. 3 calculates the velocity- 
temperature field using the relationships developed eUDOve, 

The essential features of the model as illustrated by the 
model results in Fig. 3.5(b) and (c) are that basal temperature 
decreases linearly with distance in the direction of flow (this 
assumes that heat is being conducted downwards from the fluid) . 
variation of temperature with depth involves quadratic and 
quartic terms . The depth profile of temperature has the same 
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The 
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shape at all points along the channel because T in equation 
is separable in x and y. The model has no direct ap^plication 
to the lava eSffosion problem although it does indicate what the 
temperature distribution in a fluid coming into contact with 
a conducting surface must be in order to bring convective| and 
conductive heat transport into balance. 

2 ) Simple two dimensional models - conditions for 
turbulent and fully developed flow 
In the analysis of the simple laminar flow model given 
in the last section we did not address whether the laminar flow 
model is appropriate to flow conditions experienced during 
thermal erosion and we did not consider the boundary layer 
relationships that can be expected when a flow first comes 

il 

into contact with a substrate. 

Laminar vs. Turbulent Flow 


According to Carr (1973) the mean flow velocity in a 
channel for laminar flow can be related to the hydraulic slope, 
channel dimensions, and physical pareuneters of the fluid as 
follows: o 


pgrh 


3n 


3.B.28 


where g = acceleration of gravity 
sin(j> = hydraulic slope 


r^ = hyraulic radius » CSA wetted perimeter 


n = viscosity 
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This equation is appli^cable when K (Reynolds number) < 

3 4 

1000: for common basaltic layer viscosities of 10 to 10 

poise the flow doesn't get turbulent in channels of hydraulic 
radii lOm or so until the velocity reaches about 3m/sec 
(Appendix 3.3). Solving for in equation 3.B.28. 

K /«• *m " 

^ ‘h 


» 3 X 3 X 10^ X 10^ 


3 X 371 X 10 
1 


1.000 


0.06 


6 

o 


So the steepness of slope required for the lava to flow is 

only ^out . 06®. (Carr gets 100 times this. We have not checked 

his arithmetic to see if there is an error.) In 

3 3 

his chart for n * 10 and » 10 , Carr shows the transition 

occurring at a slope of eibout 0.5° and at a velocity of about 

2 4 

3 X 10 (“3ra/sec) . However, just at n * 10 the turbulent 

transition has moved up to a higher velocity and the turbulent 

transition doesn't occur until one gets to a slope of 3.44®. 

3 4 

Thus in the range of from 10 to 10 and for hydraulic radii . 
of 10 meters or so tiie onset of turbulence is very sensitive 
to the viscosity. 

Entxx} Length - Veloc'ity and Thermal Boundary Layers 
Where a fluid first encounters a surface with a no slip 
condition a boundary layer will begin to form. The entry 

length (1_) is the distance from this point to that point 

■ 0 ' ■ 




downstream where the boundary layer thickness has grown to 
encompass the entire depth o£ the fluids 

Bntry lengths ha^;^e been determined by fluid dynamicists 
for a number of problems with simple geometry. For example^) 
for a flow into a tube the hydrodynamic entry length (1^) is 
given by; 

l^y « 0.0575 Rd 3.B29 

d is pipe diameter and Rd is Reynolds number based on that 
diameter. 

Rd - 3,B.30 

For laminar flow, Rd is <<1000 and flow becomes fully developed 
fairly close to the mouth of the tube. Analysis have also 
been conducted of flovfs over isothermal plates (Bay ley, ejt ^,11>72 ) 
which are more relevant to the geometry of the lava erosion 
problem. 

The entry lengths for thermal and velocity boundary layers 
are not identical. The Prandtl number (P) governs the re- 
lative thickness of the velocity and thermal boundary layers. 

5va P^6T 3.B.31 

The Prandtl number is solely a function of fluid properties 
and can be regarded as the ratio of kinematic viscosity v» 
to the thermal diffusivity v® =• k/(p „ Cp) 


For lava (see Appendix 3.3) the kinematic viscosity is just 
a factor of 3 smaller than the viscosity v* «* 300 Stokes. 

The thermal diffusivity of an Apollo basalt lies between 3 and 
7 X lo"^ cm^/aec* in the range 100®K to 400^K and as shown 
in Appendix 3«4 really does not increase much due to radiation 
even at high temperatures. Consequently a good value for the 
Prandtl number is: 

3 '' . u 

p » - — ^ « 2 X 10^ 

” 5 X 10"*^ 

and the ratio of the thickness of the boundary layer is ^v/^^ 
(2 X 10®)^ » 450. 


In a situation where the hydrodynamic entry length is of 
the order of 10 times the fluid depth, then the flow does not 
become fully developed thermally until It reaches a distance of some 
5000 times the fluid depth from the source. For a fluid depth 
of lOm, this distance is 50km; for 100m depth it is 500km. 
Consequently the flows carving long channels will be fully 
developed for most of their length unless they were formed 
by very deep flows. 

3) Si ^^ple two dimensional models - effects of a finite 
yield stress (Bingham properties) 

Recent research into the behavior of accreting lava flows 
conducted by Hulme (1974 , 1976) suggests that the flow of lava 
cannot simply be modelled by a viscosity. Here we review 


*Carr quotes a value 5 to 10 times higher based on the work 
of Murase and McBirney. 
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Hulme's work and its significance to the investigation of 
lava erosion* 

In an attempt to better account for lava flow dimensions/ 
Hulme (1974/ 1976) modelled lavas as Bingham fluids - fluids 
which do not flow until the stress exceeds a yield stress 
and then flow at a rate determined by the excess of the actual 
stress over the yield stress. According to this model a flow 
will rapidly spread until it thins to a critical thickness 
(Tg) at which the stress at the base of the flov/ equals the 
yield stress. Hulme showed that lava effusion rates are normally 
too high to permit significant cooling of the interior of the 
lava body before the flow comes to rest. Tensional forces in 
the more rapidly cooling skin/ even without considering its 
highly fractured character are not sufficient to restrain the 
flow. Consequently, a finite yield stress in the lava seems 
to be necessary in order to explain the finite thickness of 
flows. Even in a very thick Bingham type flow, however there 
will be heating and possibly yielding of the substrate rock 
at the base of the flow. 

If flowing lava has Bingham properties, then many of the 
relationships discussed in the last thr^f^e sections require 
modification. In addition, it is necessary to taUce into 
account an effect which is not present in the purely viscous 
theory. With a Bingham fluid the thickness of material eroded 
from the substrate must exceed the critical thickness for 
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Bingham flow (t^) befoire any nat eroaion can occur. This 
effect is illustrated schematic7;lly in Fig. 3.6. Depending on 
the flow regime, the amount of erosion taking place at th , 
base of the flow will be best represented by either the purely 
conductive relationships presented in Section A1 or the ad- 
vective relationships presented in Section A2. 

4) Simple two dimensional models - phenomenological 
model of thermo-mechanical erosion 

A plastic material is a Bingham material with a finite 
yield stress but zero viscosity. Consider a flow on a sur- 
face of varying slopes. A plastic material comes to equili- 
brium on the slope at a rate only determined by inertial forces. 

A Bingham fluid with the same critical shear stress as the 
plastic fluid, but with a finite viscosity, will take longer 
to reach equilibrium than the purely plastic fluid. However, 
it will ultimately reach an identical equilibrium distribution 
of flow thickness as the plastic fluid. To evaluate the amount 
of thermal erosion by a Bingham fluid we must be concerned 
with the finite thickness of the flow after flow has ceased 
(Pig. 3.6) as well as the eunount of material transported from 
the substrate. The plastic model can help us to understand 
the effects of yield stress on the slope profile created by 
an erosive fluid. 

In addition to the contributions to lava erosion by con- 
vective and conductive heat transport, we consider the additional 











Time <T) 


Critical thiclcness for 
Bingham flow (steep slope) 


DEPOSITION 


Critical thickness for 
Bingham flow ( shallow s lope ) 


Pig. 3.6 Thermal erosion by a Bingham fluid. Schematic 
relationships between amount of erosion and 
time for two values of the critical thickness 
for Bingham flow. Thermal diffusivity increases 
resulted in a progressive increase in slope of 
these cruves and net erosion occurring earlier 
in time. 
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mechanical effects that occ^iur when turbulence is present. 

In such situations blocks of material can be torn away and 
incorporrted in the flow. 

A simple phenomenological approach which is applicable 
when turbulence becomes important' and which includes the effects 
of both mechanical and thermal erosion is examined below. It 
depends on two ass\jmptions: 1) that there exist simple functional 
relationships between erosion rate and flow parameters. 2) 
that solution and mixing of eroded matelrial takes place on a 
time scale that is small compared to the time taken for material 
to travel the entire length of the flow. Due to rapid homo- 
genisation the . temperature of the flow can be taken to be uniform 
with depth and dependent only on downstream location; 3) that 
viscous forces (forces proportional to velocity shear) and 
inertial forces are negligible within the flow. Flow is assumed 
to be perfectly plastic and controlled by a critical shear 
stress Ty. Perfect plasticity is a special case of Bingham 
plasticity where the plastic viscosity (u ) is zero. 

First of all we will establish a relationship between 
the rate at which a flowing lava erodes material and the dynamical 
parameters of the flow. A very simple model was adopted because 
our interest is not so much to demonstrate what will actually 
occur in a given flow, but rather what is energetically possible 
for a wide range of properties of the lava-substrate interaction. 

Let us assxime a plastic fluid is doing the eroding and that 
it is a plana flow (Fig. 3.7) which has a thickness h at any 
point defined as follows: 
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pg sin© h 


3.B.33 


A i* 


where » critical basal shear stress for the plastic 

fluid 

p * density of fluid 
g » acceleration of gravity 
© * slope angle 
h » fluid thickness 

/ Iiet us also assume that the rate of erosion can be expressed 
just in terms of the basal shear stress (tb ^ t ) and the 

B p 


velocity of the fluid 

'e^B 




I 

> . C^T " i< v” 
e p 


3.B.34 


where 


^ » the rate of downcutting of the slope (see Fig. 1) 


constant 


V =* velocity 
m,n » constants 

Downcutting may be caused either by .thermal or mechanical 
erosion. 

Starting with relationships (3.B.33) and (3.B.34 ) , one 
approach to the lava erosional problem is to take an initial 
slope profile and examine how it evolves with time. Another 
approach is to solve for the slope profile which recedes but 
otherwise remains invariant with time. We will determine that 
such a profile does exist and perform the numerical integration 
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needed to define its foirin. However, we do not demonstrate 
that this slope profile is a stable profile and the asymptotic 
consequence of the erosion of any arbitrary profile. 

For uniform slope recession, the rate of downcutting 
^ can be expressed in terms of this unifosnn rate of slope 
recession (r^) 

^ = r^ sine 3.B.35 

The fractional mass gain 6f by an element of fluid thickness 
h and width Ss in the downstream dimension can be expressed as 

= (at)s “ (h|)s • v'lit'HTsT 3.B.36 

For a steady flow the mass flux will increase with s and provided 
f << 1 the mass flux (M(s)} can be expressed as: 

M(s) * ph(s)v(s) 


ph(o) v(o) */(ii 



ds 


3.B.37 


For a steady flow the plastic shear stress ( 1 ^) will also 
increase downstream due to addition of eroded debris to the 
flow. Let us assume that a functional relationship between 
X and f exists of the form 

p 

X = G(f) 3.B.38 

P 

By integrating equation (3,B.36) Xp can be'-expressed in 
terras of s. 

We shall now look at a set of flow conditions of particular 
interest when a small fractional increase in flow mass f pro- 
duces a much larger fractional increase in the plastic shear 

Stress X .In this circumstance we can neglect the increase 

"P- 











f 




in mass flux with s and therefore from equations 3.B.37 
M(s) » ph(s)v(s) ■ constant « M_ 


3. 


Substituting quantities from equations 3.B.35 and 3.B.39 in 


Equation 3.B.36 and integration: 
s > V s 

*■/(«) TTimST ^ ” f >=0 • '=/“c ^ 

o 


s 


r p/M^ f sin$ ds 


deferring to Fig. 1 we can make a change of variables from 

dH 


9 to H where sin9 

f = 


g— SO that 
s 

o . 

This simple result can be summarized as follows. Under the 


3. 


constant profile constraint for the steady flow condition the 
mass fraction of eroded material picked up by the flow can 
be simply expressed as a linear function of the elevation 
change. 

We are now in a position to solve for the shape of the 
invariant slope profile. 

3. 


a* pgh sin© 
P 


dq a C^T_^v”' 


g - resins 


3. 


3. 


TP = G(f) 


= phv 


and 


f = r. 


p/M 


H 


3. 

3. 

3. 






B.39 


B.40 


B.33 

B.34 

B.35 

B.38 

B.39 

B.40 
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We are looking for a relationship between H and sind and so 


we should attempt to eliminate the variables h, f, 

and V using these equations. 

A. Eliminate ^ between equations 3.B.34 and 3.B.35 


r_ sin© * v™ 

c 


B. 


e p 

Substitute for v in Eg. 3.B.41 

-M. 


r^ sin© 


e p 


" (;i) 


m 


C. 


Substitute for in Eg. 3.B.42 
r^ sin© = C^(pgh sinO)'^ 


D. 


Eliminate from 3.B.33 and 3.B.38 


E. 


pgh sin© = G(f) 

Substitute h from 3.B.44 in 3. B. 33 


G(f) \ n-rn 
pg sin©/ 


r^ sin© - C^(pg sin©)” (j)™ ( 

Rearranging: 

r sine . c . (G(f))'' 

c ® '• P / 


-ra 


r^ sin© = C (gM)™ • (G(f))”""^ • sin™“^© 

w © 




Sin © at 


Cg(gM) 


m 


(G (f ) 


m-n 


Finally; 


sin© = 


Ge(gM) 


m 


[(G(r^ • p/Mc • H)] ®'” 


1 

m^r 


= F(H) 




3. 


3. 


3. 


3. 










B.41 

B.42 

B.43 

B.44 

B.45 


B.46 
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In the H/ K coordinate system we can now express dH/dK 
in terms of H and Q, 


^ * tand ■ 


sinO 


(l"Sin^©)^ 


F(H) 
« * 


K . / dK . dH 


Thus K * / dK - J dH 3.B.47 

o o 

The integral in Eq. 3.B. can be evaluated once we specify 
G(f) explicitly. One form that approximates the properties 
of lava is: 

.f/fo 


G(f) * T- e 

It 


3.B.48 


Substituting in (3.B.46) 

r_ 


F (H) 


Cg(gMc) 


m 


TfO 


m-n 




m- 


3.B.49 




where 


m-n 


° Cg(gM)™ 


■po 


3.B.S0 


1 _ m-n 

Ho ~ fo '"c 


and ^ ^ (r^ • p/M) 


3.B.51 


K can then be calculated using Eq. 3.B.47. 

A prograun for integrating equation 3.B.47 (Appendix 3.5) 
allows the functional relationship between H and K to be 
determined for selected values of B, Q, and m. Some results 
appear in Figs. 3.8 and 3.9. A more advanced program is in 
preparation which allows the calculation to be made in terms 
of the more basic parameters and also calculates the other 


parameters of interest t_, f, v, and h as these vary along 

ir 


the flow. 


WsSaoLii Alt, , iaK&kts* 
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The analytical solution jjescrlbed above is only practical 
with the condition that slope remains invariant. We have 
applied numerical techniques to modelling of the evolution 
of slope profiles of arbitrary initial slope. These calcu- 
lations were performed as part of Task 4 and are described 
in our report on that task which appears in this doctunent. 








Horizontal Distance (K, 



Horizontal Distance (K, km) 



Fig. 3.9 Lava flow t>ed profiles for steady state retreat of a lava 
flow eroding a substrate. 
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TASK 4 LAVA EROSION: LABORATORY 

SIMULATION AND TWO DIMENSIONAL 
FINITE ELEMENT MODELS 

The two aspects of this task were conducted in paralloi 
beginning in June 1979. The laboratory wax model experiments 
were very limited in scope and were concluded by the middle 
of August 1979. Development of the two dimensional finite 
element codes for modelling thermal lava erosion continued 
up to the end of the contract period. 

A. Wax Model Experiments 

Beginning in June 1979 we performed a series of experi- 
ments simulating lava .erosion by physical modelling methods. 

Prof. Jay Melosh, at that time of Caltech, permitted us to 
use laboratory facilities set up for conducting physical 
simulations of ocean floor spreading using hot wax. We al- 
located some discretionary funding for -equipment to investigate the 
feasibility of simulating lava erosion using hot wax. Mr. Jose Helu, 
Prof. Melosh's former research assistant in the Caltech physical 
modelling experiments, assisted us with these wax model ex- 
periments, The results of these experiments exceeded our ex- 
pectations and proved to be a valuable guide to our theoretical 
modelling efforts. 

1. Experimental Approach 

To simulate erosion of rock by a flow of molten lava we 
pumped hot wax into a channel of rectangular cross sections 
which had been formed in a wax substrate.. The experimental 
set up is shown in Fig. 4.1 (a) (b) . It includes a thermomix 
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heater to heat the wax to a uniform temperature, a pump to 
supply the fluid at a uniform rate, a wax casting with a channel 
of rectangular cross section, and a pump for collecting wax 
reaching the end of the channel. Two experiments have been 
conduqted so far using poly thy lene glycol (Carbowax 4000) for 
the fluid and substrate. 

2. Results 

Photographs of the wax casting after one experiment are 
shown in Fig. 4.1(c). In this experiment, material 

was eroded along the entire length of the channel. At the 
inlet the hot wax dropped onto the channel floor and this 
resulted in the formation of a depression; at the outlet hot 
wax cascaded over the edge of the channel and eroded headward. 

An interesting feature of the erosion process is evident 
in the channel cross sections (Fig. 4.1(d)). There has been 
significant removal of material by lateral erosion excep-^. 
near the surface of the flow where hot wax has solidified and 
accreted creating a significant overhang. This behavior was 
not suggested directly by the theoretical modelling efforts 
of Carr (1974). However, the forms of the isotherms developed 
when a hot fluid occupies a cavity of rectangular cross section 
(Carr, 1974) suggest that thermal erosion will develop more 
rapidly below the surface of a flow than right at the surface, 
so in retrospect, the results of the hot wax erosion experi- 
ment are not surprising. 
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The depth of the channel before and after the experiment 
was measured with a digitizer and the total erosion at the 
channel center was calculated and plotted for one of the ex- 
perimental runs (Fig. 4.2). The large amount 6f erosion 
near the channel source was a consequence of a 'waterfall' 
effect created when the hot wax first spurted on to the sur- 
face of the channel. Between 10cm and 35cm the aunount of 

>'[ 

erosion declines steadily as the hot wax cooling through 
transfer of heat to its environment and assimulation of ori- 
ignal coq,l wax from the chamnel wall and floor. Channel 
sections were also digitized illustrating the variation of 
erosion/and accretion as a function of depth and are displayed 
in Fig. 4.3)., 

In Table 4.1 we have summarised the effects of thermal 
erosion by the hot fluid wax on the wax substrate. A, signi- 
ficant parameter is the transport and erosion efficiency, 
i.e., the proportion of material that flows out of the channel 
that is derived by erosion of the’ substrate. It is most con<- 
veniently expressed as the ratio of the volume of material 
eroded to the volxime of material which flowed in the channel. 
Material removed due to headward retreat at the channel out- 
let was excluded from this estimate. V 

The relevance of hot wax simulations to the behavior 
of lava flows has been considered by Hodgson (1969) (Appendix 
4.1) and at a workshop held at the Los Alamos Scientific Lab- 
oratory (Widdicombe ^nd McGetChin, 1976) . The emphasis was 
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DATA ON THERMAL EROSION EXPERIMENTS WITH HOT WAX 



on the emplacement and mechanics of a flow of lava and not 
on the erosion of a substrate (substrates which did not melt 
were used ) , however/ the similitude arguments developed (see 
Appendix 4.1) are still pertinent here. Thus, the approach of 
using scaling relationships to describe the results in lab- 
oratory experiments should also work with the erosion problem. 
The wax we used (CarbowaJC 4000) was the saune that was used in 
experiments for the lava experiment at Los Alamos. It melts 
to a clear liquid at 55°C, forming a non-Newtonian fluid. 

Many of the similitude numbers match those for lava (see Ap- 
pendix 4.1); moreover it does form analogs of lava flow features 
such as natural levees, lava tubes, and pahoehoe surface 
textures . 

3. Further Investigations 

Laboratory experiments with hot wax provide a powerful 
tool for developing qualitative and quantitative understanding 
of the mechanism of thermal erosion at a geological scale. 

Such experiments can be used to validate numerical methods 
for calculating lava erosion; study thermomechanical erosion 
in circumstances which are not practical with numerical methods , 
e.g. , pulsed flows, channels of varying cross sections, 
sinuous channels, and study the formation of morphological 
and structural features of the lava erosion regime. 









Specific goals of future investigations with hot wax \\ 
ntodels might be: 

1) To validate ^e numerical codes for thermal erosion 
developed for slah flow and channel flow as de-> 
scribed in Task 4B using laboratory experiments 
with flows of hot wax. 

2) To examine the transition from erosional to de- 
positions! conditions using the physical modelling 
approach . 

3) To investigate how the efficiency of erosion de- 
pends on such parameters as temperature, viscosity 
and flow rate. 

4) TO examine the effect of pulsed or intermittent 
flow on erosion rates in a channel. 

5) To examine channel cross section and profile de- 
velopment under a variety of flow conditions. 

6) To examine flow in sinuous chcuinels. 

7) To apply the validated theoretical codes to the 
modelling of large scale thermomechanical erosion 
on planetary surfaces and to investigations of the 
origin of certain classes of erosional features on 
planetary surfaces. 

An improved experimental set up would be needed to achieve 
the accuracy and reproducibility needed for meeting these 
objectives. Equipment additions would be required for: 


1) improving control of the flow rate and temperature 
Of the fluid wax entering the wax channel, 

2) improving monitoring of the temperature of the wax 
substrate, and of the fluid wax along the length 
of the channel, 

3) refining measurement of channel profile and cross 
section, 

4) making additional molds to permit more experiments 
to be run. At present, each mold takes several 
hours to coool and this limits the number of ex- 
periments, 

5) use of channels with varying amounts of sinuosity 
in plan, 

6) using channels with obstructions and constrictions 
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B. Two-Dimensional Finite /Element Models 

Our approach to the problem o£ modelling thermal 
erosion using finite element techniques involved two 
stages. First, we spent some time developing a basic under- 
standing of finite element techniques and the stability 
relationships that restrict their applicability. Then, we 
embarked upon a serious formulation of the lava erosional 
problem and developed a series of computer programs of 
progressively increasing complexity. 

1. Modelling Methods 

Our initial investigations were designed to gain 
familiarity with some commonly used finite difference 
algorithms. We exercised these algorithms on some simple 
thermal problems and compared the observed convergence 
properties with those anticipated theoretically. 

The vorticity-stream function method was anticipated 
being used for our numerical investigations of thermal 
erosion. In this method the momentum equations of Navier 
Stokes for incompressible flow are cast into the form of 
a parabolic vorticity transport equation. This equation must 
be solved for tiie specified boundary conditions and the 
solution used as the inhomogeneous term in a Poisson 
equation for the stream function. 


After solving this elliptic problem for the stream 
function by iterative methods, velocities are obtained as 
the partial derivatives of the stream function with respect 
to the spati/jil coordinates. The vorticity transport 
equations are identical in form with the heat transport 
equations and so this approach appears to be attractive 
for studying thermal erosion where heat and vorticity 
transport are coupled. In a problem where there is no flow, 
or where the flow is specified, the solution of the thermal 
analog of the advective diffusion equation for vorticity 
transport represents a complete solution to the problem of 
interest. 

Several finite difference schemes were examined for 

solution of these equations. The Dufort-Frankel leapfrog 

is a single step, three level finite difference method for 

solving the advective diffusion equation. This method is 

unconditionally stable as there is no limitation on the 

value of the diffusion number (d * This means 

A>C 

that the mesh size can be reduced without requiring a much 

larger reduction in the size of the time step. The only 

steibility restriction is the invisCid requirement on the 

courant number (C = . of C < 1. We developed a 

Ax 

computer program to apply the Dufort-Frankel method to a 
two-dimensional thermal diffusion problem as described in 
Appendix 4.2. 


Wllbh the experience gained in these exercises v/e were in 
a position to address the specifics of thermal erosion. In this 
we benefited from the guidance of Stephen Keihm. Our approach 
here was to investigate models of progressively increasing 
complexity. The models are described here below; the computer 
programs to implement the models are included in Appendix 4.3. 

2, Thermal erosion for a plane laminar constant viS" 
cositv flow with a yield temperature. 

This model resembles the approach originally formulated 
in a one dimensional form by Carr (1974) in that when the 
substrate temperature reaches or exceeds the yield temperature 
the material participates in the flow. 

Approach 

We first assume the velocity field for a plane , laminar , 
constant viscosity flow: 

V(y) » ^ sinij) (H^ - y^) 4.B. 1 

where H *= flow layer thickness 
\> » kinemetric viscosity 
g * gravity accel. 

Then, as a first approximation, with the velocity field 
given , the problem can be treated from a purely thermal approach. 
We have used a Forward Time Centered Space (PTCS) algorithm 
to model the conductive heat flow; heat is advected at the 
fluid velocity using a Forward Time Forward Space (FTPS) 
algorithm. 


Pig. 4.4 Geometry for Finite Element Simulation of Thermal 
Erosion. 





For the rectangular mesh as shown abovr 

- T? . 

pCAVAy • pc ' - V^^jAypC - 


• Ay( ^i“l>j " + kAy 

' AX ' ' Ax ' ' 


tie ^ AX ^ + kAy ( ij.) 


4.B.2 


But V. , j * V. . since V is indepent of longitudinal direction 


- T?” V 

Aij-Z— If-i * -i ^'^i-l,j " "^irj^ 

At Ax 


4.B.3 


+ <Vi,j • ^'*i,j * '^i+i,j> 




Where * V(y) is given adsove. 

Since the velocity field is assumed fixed, the removal 
of material from the supporting wax layer must he evaluated 
by a thermal criterion. Two approaches might be used on 
limiting cases: 
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1. Evaluate the depth of erosion as a function of the 
downslope coordinate / x, by fitting the temperature 
profile of tlie wax base to a quadratic function. 

The depth of erosion would then be estimated by the 
depth to which a 'yielding' temperature, Ty, was 
attained. This is a similar approach to the one 
dimensional model with conduction and convection 
described in Section 3.B.2. This would correspond 
to a lower limit estimate of the erosion since it 
does not take into account the gradual removal of 
material and resultant new 'contact' of underlying 
base material to tlie high, flowing fluid tempera- 
ture. 

2. Add a stipulation in the numerical model that as 

soon as any base reaches the yield temperature, 

Ty, it is replaced by the adjacent fluid temperature. 
This is a similar approach to the 'modified' one 
dimensional model with conduction and convection 
described in Section 3.B.2, It is also similar to 
the scheme used by Carr C1974) and should represent 
an upper limit on the erosion with heat of melting 
effects neglected. The aunount of erosion is again 
calculated by the depth of the yield temperature 

in the wax base. 
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Boundary Conditions 

1. Boundary nodes of the wax base radiate to ambient 
space as do the surface modes of the wax flow. 

2. At the source of the wax flow (x » 0; 0 £ y £ H) 
an input ('left to right*') flux corresponding to 
the hot wax temperature is used. 

3. At the 'spill end' (x » L; 0 £ y ^ H) , zero con- 
ductive flux is assumed. 

Since the flow velocity is assumed constant, we start 
with all nodes for which y i H set to the hot wax temperature, 
with the wax base ambient throughout. 

There is an input at the left of T^, and the output 
boundary is held at conditions of no heat flow into or out 
of the boundary, as is the bottom of the solid wax. The top 

4 

surface radiates at , but again, no convective or advective 
heat flow into tlie boundary is permitted. These conditions 
are realized by creating a dummy set of values outside the 
boundary which are set equal to the values on the inside, 

4 

Causing no heat flow across the boundary. The condition 
is dealt with in a separate step. This leaves only the solid- 
liquid interface to be dealt with. 

In earlier test versions , the interface position was held 
fixed — -i.e. irrespective of the temperature, liquid never 
solidified, and solid never melted. This is somewhat equi- 
valent to having a thin insulating sheet at the interface, 
and though unrealistic, served to test the algorithm. 






The present version allows for melting and solidification 


in the following manner: 

a) A threshold temperature is defined at which 

fluid flow will occur. 

b) Solidification is assumed to begin at the far end, 
and erosion is assumed to begin at the near end. 

c) If a cell which was in the original fluid regime 
drops below the threshold, temperature advection is 
turned off. Further, the heat that would have been 
advected into the cell (assuming, of course that the 
cell in front of it is still fluid) is advected 
straight up into the cell directly above it. 

d) If a cell which was in the original solid regime 
rises above the threshold temperature, the fluid 
is 'mixed' with the fluid directly above it. This 
is accomplished by advecting a quantity of heat less 
than or equal to the amount of heat needed to bring 
both cells to the same average temperature. It is 
presently set up to mix completely in one second 
(100 iterations). 

e) The extra heats needed to accomplish this are held 
in a separate array and added to the cells (or sub- 
tracted from them) at the next iteration. 

The FORTRAN program for implementing this simulation is 
reproduced in Appendix 







52 


Resul fcs 

■Jhree runs of the program were performed with identical 
parameters excepted for velocity. %ese runs simulated the 
effects of a relatively slow moving lava flow (V = lOcm/sec) 
V.m^ fact flow (V ^ lOOcm/sec) . The parameters used 

in these runs and the number of iterations performed and plots 
generated are shown in Table 4.2. 

The immediate Output of the program was a listing of 
temperature values on a 12x17 grid of points (eg. Fig. 4.5(a)). 
More than 200 printouts of this format were generated. In 
order to better visualize the results of the experiment a 
graphical representation of the same data as isotherms was 
generated (eg. Pig. 4.5(b)). by interpolating individual 
temperature values from the printouts. Only a small number 
of these plots were generated. 

In Pig. 4.6 we show the isotherms for each run at a series 
of times into the experiment. The horizontal dashed line in 
each plot represents the original interface between the sub- 
strate at 20°C and the fluid at 1300°C . The 600°C isotherm 
which is the lowest one plotted provides a useful reference 
to indicate both the flow of heat into the substrate and the 
progress of erosion of the substrate. 
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In Run 1, the 600®C isotherm progresses by about one cell 
into the substrate during the 96.7hr experiment except close 
to the source where progress is much more rapid. In Run 2, 
the 600°C isotherm advanced by all sizes into the substrate 
during the 96.1 hr experiment at the downstream end and by 
larger amounts nearer the source. In ’contrast with Run 1, the 
isothem maintains a significant slope for the entire length 
and is not as sharply curved near the source. In Run 3 the 
600°C isotherm advanced by 3h cell sizes into the substrate 
during the 96. Ihr experiment at the downstream end. 


TABLE 4,2 PARAMETERS FOR TWO DIMENSIONAL 

LAVA EROSION SIMULATION 

KAPPA * 0.04 
ROC * 0.25 
H = 0.50 

LENGTH » 1.0 X 10^ 


G » 980 
ANG * 30 





iTfwrjoN* 

1 2 

1 t)«9.0n 
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Fig. 4, 5 (a) Matrix of temperature values 
for Run 3 (VMAX * 10 00 cm/sec) 
at 96.7 Hrs. 
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Fig. 4.5(b) Isotherms for data shown in Fig, 4.5(a). 
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II 

CONCLUDING REMARKS 

Material from this report is part of a paper which 
is currently in preparation on the geology of martian 
rille systems. We envisage that future work on this problem 
would logically involve a joint experimental/theoretical 
attack on lava erosion. A proposed follow-on activity 
submitted in July 1979 (Appendix 4.4) which would have under- 
taken such a project was not recomended for funding and at 
this time we have no plans for continuing investigations 
of this problem. 
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MODELS OF EBOSION RELEVANT TO 

QIANNEL FORMATION ON MARS 


by 


James A. Cutts 


Planetary Science Institute 
283 S . Lake Ave . , Suite 218 
Pasadena, California 91101 


ABSTRACT 


Evidence for the role of lava erosion in the formation of some of the 
martian channel features between western Chryse Planitia and Lunae Planum 
is presented. EarlJLer numerical modelling studies of thermal erosion by 
lava are reviewed and deficiencies in the model are identified. An alternati.ve 
simple one-'-dimensional model is presented which rectifies some of the 
weaknesses in the earlier model. However, this model also fails to ade** 
quately characterize the lava erosional process, ihe features of more 
sophisticated numerical models which are needed to fully depict thermal 
erosion by lava are described. The relevance of these studies to the 
turbulent regime which may apply in real«>world lava erosion are assessed. 


For presentation at the 

pla;ietary geology field conference on basaltic volcanism 

Snake River Plain, Idaho 
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IMTHODUCTION 


Imges obtiained by the Viking Orbiter have reveaied that many of the 
so~called naxtian chai'inels occut in an intimate association witti plains 
vulcanism. Ihese discoveries contrast with the impression provided by the 
earlier Mariner 9 orbital Imgery that the channels were older than the 
pl 2 uLns units. We have studied a group of such channels in some detail 
(Pig. 1) . The results of this investigation are reported in an article 
recently submitted to Icarw (Cutts and Blasius, 1977a) . The conclusions 
of that paper are that these channels have been formed by erosion by the 
overflow of lava from Lunae Flwum into Chryse Planitia across the 
intervening mountainous divide. In this short note, I wish to merely 
summarize the key points of the evidence bearing on this interpretation 
and to- then report on some recent theoretical examination of the process of 
lava erosion. Field studies associated with this conference will add some 
further insights into the mechanism. 

EVIDENCE FOR lAVA EROSION BETWEEN LUNflE PLANUM AMD WESTERN CHRySE PLANITIA 

Predicting precisely what a large lava erosional channel would look 
like involves some detailed morphological arguments which are too con^lex 
to repeat here. There are two key pieces of information however, which 
support lava erosion origin for the channels. 

1. The cheuinels connect two surfaces which are identified as volcanic 
plains by conventional photogeologic criteria. 

2. ©le ages of these channels determined by two largely independent 
methods ~ small crater density and large crater superposition/ 
intersection relationships - are indistinguishable from the ages of 
the adjacent plains. The first method yields an age differential 







of le$i than 0.22 at the 95% confidence level; the second 

method an age differential of less tl^ 0.24 at the 95% confidence 

level. 

Penec<^ntenporaneity in age supports a volceuiic origin for the rilles 
through two different lines of argument. Firstly, it is a necessary condition 
for both pledns and rilles to have been formed by lava. Secondly, it requires 
a remarkable coincidence for an erosional episode involving some other high 
density flowing medium to have occurred just once in M£u:s history and almost 
contemporaneously with the effusion of plains basalts. 

lAVR EROSION MECH .fl.N ISMS 

The importance of lava erosion as a possible agent of valley foisnation 
Cn lunar and planetary surfaces has been stressed by Carr (1974) . Greeley 
and Hyde (1972) have described field evidence for lava erosion on the earth. 
Carr (1974) has attenpted to model the process of lava erosion in a laminar 
flow of lava« 

Stimulated by Carr's investigations we have attenpted to understand the 
mechanism of lava erosion in more detail. ' We have re-examined his model and 
found certain deficiencies. We have performed some preliminary calculations 
of our own in an attempt to improve the model. However, we have concluded 
that the problem is such that only a full three dimensional treatment C 2 in 
accurately protray erosional phenomena in a lava chamnel although 
significant insights can be gained into the laVa erosion process by two- 
dimensional models of a planar flow in which one dimension is along the axis 
of flow. Neither OVrr’s model nor our first-order model incorporate this 
dimension. 

Carr's model of the mechanism of lava erosion 

The only detailed investigation of a mechanism of lava erosion that 


has been repainted is the work of Cairr (1974) . Mechanical and thermal 
effects both occvur in lava erosion, but only thermal erosion is amenable 
to mathematical modelling. At the base of a lava ^ yhere the downcutting 
process takes place, thermal effects will dominate mechanical effects, 
since the floors of lava tubes are smooth and glazed providing li,ttle 
opporttinity for mechanical plucking. 

Carr (1974) has considered thermal erosion as a three step process: 
first, the wall is heated, then it flows, finally it becomes incorporated 
in the lava stream. Carr asserts that the first stage operates independently 
of whether the flow is turbulent or leuninar. Ihe second stage, flowage 
of the walls , depends on viscosity and tangential stress . The viscosity 
depends only on the temperature in the walls; the tangential stress only on 
stream depth fob both laminar and turbulent flow. I£ turbulence is a 
factor , then it controls incorporation of the flowing wall materials into 
the stream. 

The lava erosive mechanism has been modelled by Carr in order to 
estimate erosion rates. He assxomes that a lava channel grows by heating 
of its wall materials until they reach a yield tenperature (Ty) at which time 
they become sufficiently fluid to flow and become part of the lava stream. He 
determines the temperature in the vicinity of a lava channel by using the 
two dimensional heat equation which he solves by applying a stcuidard relaxation 
technique to a variably sized array of points representing temperatures in 
a cross section through the channel and its surroundings. As a starting 
point, the dimensions of the channel are defined, all points within the 
channel are set at the lava teirperature and the rest of the points are set 
at 0°C- Two further conditions are that the surface except for the channel 
is kept at 0°C 2 Uid that when any tenperature reaches a yield teiroerature it 
is replaced by the lava tenperature. This second condition inplies that 
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wall material reaching the temperature ty becomes incorporated into the lava 
stream. 

Calculations made with this model suggest that thermal erosion rates 
are very sensitive to the .difference between the lava temperature and the 
yield temperatiire. For temperature differences between 0°C and 100°C the 
erosion rates vary from zero to approximately one meter per month. However, 
increasing the ten^ierature difference by a further 200°C does no more than 
quadruple the erosion rate to approximately one meter per week (Fig. 2 ) 

These results, although certeiinly within the reange of erosion rates 
that can be observed in nature, may contain serious errors because of 
deficiencies in the model. !B!iese deficiencies occur in three areas: in 

the representation of the temperature-velocity field at the base of the 
flow; in the characterization of the temperature and stress dependent 
viscous properties; and in the algorithm vdiich models the thermal erosion 
process. deficiencies are a consequence of achieving mathematical and 

conputational tractability with a 2-D model of channel flow. 

The first two problems with the model are coupled. Rather than solve 
the Dodifie’d Navier-Stokes equations for velocity amd temperature, this 
conrolex problem of interacting velocity and temperature fields was sinplified 
by separating the transport of mass from the transport of heat. To effect 
this simplification it was necessary to aissume that: the material in the 
substrate of the flow abruptly undergoes a transition from a rigid solid to 
a mobile fluid at a yield temperature Ty. However, laboratory and f.teld 
!neasuxeIl^ents on molten lava as Carr discusses, reveal a gradual, rather than 
an abrupt change from a rigid solid to a mobile fluid over a temperature 
rarige of up to several hundred degrees. There may be materials or circum- 
stances for which this yield temperature is a valid approximation but this 
has not been demonstrated so we believe some skepticism is justified. 
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An even more serious difficulty arises withi the , algorithm that is used 
to calculate erosion rates from the propagation time of a temperature waw 
into the channel wall. While this algorithm is only applicable to the 
yield- temperature formulation of lava properties it appears to give misleading 
results even with substances for which this yield temperature model is a 
' valid approximation. As we shall now show^ the rates of erosion predicted 
are critically dependent upon the characteristic dimensions of the array 
of points used in the numerical solution for the propagation of 
temperatiire into tb.e chaiuiel wall. 

Let the eurray of points upon which temperatures are represented have 
a characteristic separation at the channel wall when the i^ cycle 
of iterations. Which results in an erosional event, takes place. time ?ti. 

for the temperature at the first array point inside the channel wall to 
reach the yield tenqaerature , thereby bringing about this i^*' erosional event, 
will scale as the square of the array separation. 

,^'ti = Ki ^xi^ 

where Ki is some function of the thermophysical properties of the medium and 
the thermal conditions at the onset of this ith cycle. 

The total time tjj for the yield temperature to propagate successively past n 

array points is an arithmetic sum of the tx, 

n 

i*l 

and it will propagate a total distance Xr J^spresenting the total amomt of 

wall material eroded given by 

n 

i»l 
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TO find the behavior of tn as the array separations bocoroe infinitesimally 
small, we can assune that all of the ^xi an<?, all of the Ki are of the same 
order as one another. Then substituting in equations (2) and (3) we find that 

t Ki • 5 xi * 5«n (4) 

n ^ ^ 


Thus, the tuns to erode a given thickness of channel floor depends on the 
characteristic array separation which is a mathematical artifact designed 
to aid in the solution. Moreover, as the array separation tends to zero#, 
which should lead to inproved precision in the solution, the time to erode the 
same sickness of chaumel floor tends to zero also. It appears that erosion 
times calculated by .the model may therefore be largely an artifact of the 


particular grid separations chosen. 

We consider that to achieve realistic estimates of the rates of lava 
erosion it is necessary to explicitly consider variations in temperature 
and velocity ^ the flow direction as well as perpendicular to it. A.s a 

full 3-dinensional treatment becomes extremely elaborate and costly, an 

» 

alternative is the study of planar rather than channel flows- This will 
illuminate the basic physics at the base of the flow and will also provide 
useful, better than order of magnitude, estimates of erosion rates given 
that effects of channel walls must be ignored. It is probably that estimates 
of parameter interrelationships will be accurate to a factor of 2. 

Not only would a 2-D planar flow treatment yield valid erosion rates 
but it also may yield numerical results that can be compared with dimensional 
information on liinar features formed by lava erosion such as topographic 
profiles along lunar rilles eind, excavated volumes of lunar rilles. The 
important transitional region between lava erosion ^uld accretion could also 






be studied. However « some problems such eus channel profile and sinuosity 
would obviously be beyond the range of applicability of the plauiar flow 
treatment. 


An alternative model of lava erosion (Model 2) 

We have made some exploratory calculations of temperature amd velocity 
fields within planar flows. In one set of calculations, detailed below, 
we studied the flow in the substrate beneath a planar lava flow assuming 
the relationship between viscosity and temperature for lunar rocJcs given 
by Cukierman (1973) . It represents anotlier approach to the erosional 
problem by Carr (1974) but we en^hasize that it suffers the same fundamental 
problem, of ignoring the downstream dimension. We adopted atn essentially . 
arbitrary criterion for the depth of erosion, namely that depth at which 
the flow velocity exceeds 10 cm/sec. Our initial objective was to determine 
the rate of erl'jsion, i.e., the variation of this depth of erosion with time. 

It proves fruitful to consider separately the dynamic behavior of the lava 
flow and that of the rock substrate. In the only numerical simulation 
completed to data; we have considered a steady state flow of an isothermal 
lava with infinite heat capacity and conductivity over an inclined rock 
substrate, focusing on the behavior of the rock substrate. Temperat\ire emd 
viscosity versus time profiles for the rock substrate were developed. 

Prom these profiles, the rate at which melting rock joins the flow, and 
hence how the depth of melting increases with time have been determined. 

The motion of the substrate beneath a planar lava flow was modelled adopting 
a simple functional relationship between the viscosity of the substrate and 
its temperature. We assume that the lava channel grows by heating of the 
substrate such that the substrate mobilizes ai'id begins to flow in the direction 
of the lava stream. This flow is assumed to be laminar and the substrate 
is assumed not to mix turbulently with the free lava stream. However, in 
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order to aimplify the problem we must also tidce the lava itlow to be an iso- 
thermal reservoir of heat at temperature M we are concerned only with 

the behavior at the base of the flow, a one dimensional solution for a 
planar rather than a channel flow is adequate to expose the basic physics 
of what happens in a channel flow (Fi?* 2) . 

As a starting point we ass'jme that the" tenperature of the substrate 
is To and the thermal diffusivity is K- At time t=0 an insulating layer 

separating the lava from the substrate is removed bringing the two into 

• * 

thermal contact. The tenperature within the subs.^ate has a straightforward 
analytic ^‘<* Aion if K is independent of temperature. Thermal diffui^vities of 
silicate mH<:irials in fact, show considerable variation (a factor of 10) 
over the temperature ranges in question but we will ignore this coitplexity 
in this treatment. Thus: 

Ts - (Tl - - erf ( I /"T . 's’ 


The 


viscosities of lunar basalts vary from i>erfact rigidity on 


time scales relevant to the duration of a lava flow^near below 900 C to 


approaching the consistency of motor oil, at 1200 C. We have performed a 
parameterization of data by Cukierman et , (1973) for the viscosity(y of a 
basalt which is considered here to form the bed or substrate of a chauinel 


flow. 


( 6 ) 


log = A + B/T + C/T 

7 


vhere A >= 3.145 


B = -1. 84 X 10 
.7 


37 x 10 






c 


The tesperature within the substrate will vary sufficiently slowly for inertial 
effects to be neglected. Consequently the equation of inocion for the 

i| 

fluid JLs : 


dy 






(7) 

»♦ 


Integrating , 

T'- + yfis> 9 sinG 



where y =* 
u(y) = 

^ (y) = 

h “ 

Pl “ 


depidx beiow lava/substrate interface y 

velocity of substrate parallel to Interface 

viscosity 

depth of lava 

denSJ.ty of lava 


Pg » density of substrate 
g * acceleration of gravity 
0 » slope magnitude 


( 8 ) 


\ 


A solution consistent with the boundary equation can be obtained by integrating 


( 8 ) 








g sin0dy 


(9) 


The rate of flow along the channel is estimated by integrating (9) 


dy 


( 10 ) 


Numerical solutions for the flow rate as a function of time have been 


obtained and are shown for a variety of lava temperatures in Fig.H* 


The 
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figure clearly illustrates how the substrate flow builds with time. Bie 
semi'- in finite character of the boundary conditions results in no net removal 
of substrate material. For an estimate of the deepening of the channel 
produced when all the lava flow is drained we have somewhat artificially 
taken the 10 cm/sec velocity level as the base of the flow, figure S 
illustrates how this level changes with time. 

Our results, identified as Model 2, are compared with the thermal erosion 
data obtained by Carr in Figure 2. Also shown is the '/iscosity- temperature 
reladonshlp used in our calculations. Model 2 must overestimate the time 
required to remove one meter of substrate eus it Incorporates no means of 
bringing fresh hot lava into contact with the substrate. Howe'ver, we 
en^hasize that we consider neither erosion rate vs. temperature relationship 
in Fig. 2 is valid because of the inadequacies of both models. 


DISCUSSION 

Although the calculations made above illustrate how the viscous properties 
of the substrate rocks can be brought into the assessment of flow processes 
occurring during lava erosion, they do not permit a meaningful calculation 
of erosion rates and they provide no insights into the changing chcuracter 
of the flow regime along the channel. In order to do this it is necessary 
to set up a 2-dimensional model treating the temperature and velocity fields 
in both the lava flow and the substrate and taking convective as well as 
conductive heat trzuisfer into account. 

The goals of such an investigation would be fourfold: 

1. To determine rates of lava erosion with realistic numerical models 


of heat and mass transfer processes in the bed of a planar lava flow. 

2. To estimate the efficiency of lava erosion (the ratio of eroded 
material to flow material) as a function of slope, viscosity and substrate 
material properties. 


pms If 
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i. To oxamine the temporal evolution of the profile of the bed of a 
planar lava flow due to lava erosion as a function of the initial profile,, 
the properties of the eroding flow and the properties of the substrate. 

4. TO study the transition between lava erosion and lava accretion 
as a function of slope profile, the properties of the lava fl'ow, and the 
properties of the substrate. 

Hxe geometry of the model for simulating lava erosion of an inclined 
substrate by a planaur flow is depicted in Fig. 3. Heat and mass transfer 
between and within the flow and the substrate would be modelled using 
the general equation of heat transfer (Landau and Lifshitz, 1959, p. 185) 
modified to handle temperature-dependent viscosity (Bay ley, et al., 1975). 
Thexindynamic and kinetic variables would be computed on an array of points 
within the flow materials and within the substrate. 

Initially, the substrate is assumed to have a uniform temperature 
Tso and the lava flow a uniform temperature they are assumed not to 

have been previously In thermal contact. Initial velocities are zero in 
the substrate, uniform and non-zero in the lava at a rate appropriate to 
the constamt viscosity and the slope. 

At t==0, heat excahnged is initiated across the half-plame downstreaun 
of a central reference point. Ihis artifice permits investigation of the 
profile of lava erosion pecullair to the origin of the flow ais well as the 
more gradually chainging gradients for downstream. Initially only conduction 
is significant but as the thermocynamic structure of the flow field evolves, 
convective trauisfer will become more significant in heating they substrate. 

Ixperimentally determined data on lunar rocks would be used for the 
thernDphysicail properties of the lava flows amd sixbstrate. Some of these 
paraz£ters such US viscosity and chermal diffusivity are temperature 
dependent. To more precisely express the phenomena of flow in geologic 
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aatesials^ a fom of plastic flow model is needed such as the Bingham fluid, 
in which the deviatorio strai.n rate is proportioned to the deviatorio stresj 
difference in excess of a chosen yield stress. 

How relevant are investigations such as these when real lava erosion 
problems Involve channel flow not planar flow and turbulence is likely to 
play an in^ortant part in the erosion process? 1!he same response suffices 
to answer both questions, ihe ^action* in thermal lava erosion^ as far as 
the basic issues of incising channels into a cold rocky substrate is concerned^ 
takes place at the base of the flow. For modelling this interaction in a 
channel the planar flow model is quite adequate for providing initial insights 
even though it requires modification for accurate estimation. At the 
base of the flow also^ turbulent effects are at ''t minimum. Certainly it is 
necessary to understand the role of turbulence in channel formation as it 
certainly has an ispact. But we first of all need to understand the processes 
that occur in less complex flow regimes . 

SUMMARy 

Possible lava erosional channels on Mars have been described. Existing 
work on the modelling of lava erosion has been reviewed and future 
directions for numerical modelling studies of lava erosion have been identified. 
Field studies of prehistoric lava erosional features such as those present in 
the Snedce River plain will provide important insights into the process of 
lava erosion on lunar and planetary surfaces. 
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FIGURE CAPTIONS 

1. Martiem channels Vedra Valles /Maumee Valles (top) cmd Maja Vallls (bottom) 
Vedra, Maumee and Maja Valles are interpreted to be lava erosional 
features. Maja Vallis also exhibits later modification by large scale 
flow of some substance whose nature is presently \£|’Jcnown. 

2. One dimensional model for planar flow of a viscous substrate beneath an 
isothermal lava reservoir. 

3. Estimate of time to erode one meter of substrate with lava as a function 
of lava temperature. Model 1 is due to Carr (1974) . Model 2 is described 
in text. The third curve is the viscosity temperatiire relationship 
for a lunar basalt. 

4. Volume rate of flow of substrate j)er unit width of a planar flow as a 
function of time (t) and lava temperature (Tj,) The increase with time 
represents both a deepening of the depth of substrate that flows and 
an increase in velocity. 

5. Estimated depth of erosion as a function of time. ‘for Model 2. 
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TEMPERATURE DEPENDENCE OP 
ROCK CONDUCTIVITY 

The thermal diffusion calculations presented in Appendix 
3.1 do not incorporate the temperature dependence of rock 
conductivity. Here we describe this dependence. 
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THERMAL CONDUCTIVITY OF ROCKS 

1. Ambient temperature values 

Paper by Harai and Winkler (1974) lists the thermal diffusivity of an 
Apollo basalt (porosity 10%) as 

7 X lo”^ cmVsec at 100°K 
3 X 10 ^ cm^/sec at 400°K 

For purpose of checking units we will assume an average of these values 
K 5 X lo”^ cmVsec 

Thermal diffusivity may be converted to thermal conductivity using density and 
specific heat so that thermal cond of lunar basalt (kj^g) is given by: 

lb -3 • 

k2ie - K (9 c ; P 5 3 g cm 

C 0.2 cal /g/°K 

» 3 X 10 ^ cal cm ^ sec ^ ^ (1) 

This can be compared with the value for granite on PE-16 of Hemdbook of Physics 
and CSiemistry which 

kgj. = 1 to 3 k cal m ^ hr ^ deg ^ (2) 

. -3 

Conversion factor to units of (1) is 2.77 x 10 

Therefore thermal conductivity of granite is 

K »3to9xl0^ cal cm ^ sec ^ ^ 

9*^ 

2. Dependence of conductivity on temperature 

The conductivity of silicates shows a significant temperature dependence 

including cube law and inverse terms. The calculations on the next page use 

the relationship due to 

k = A + B , + CT^ 

'T 

, _3 2 

where A = 015235 X 10 cm /sec 

B = 0.1561 cm K/sec 

C = 0.1367 X lO'^^ cmVsec 
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Over the temperature’ range 100°K to 2000°K however^ the actual variation in 
conductivity v/ith temeprature is not large as the B and C terms are not large. 
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Elailfc^wavis vclocttfcs and ihcrmni df/Tusfvirks of Apollo 17 rocks 
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Lunar Anorthosite 
77017, 24 
(Ts 400*K) 


10 ^ 10 10 *' 10 *^, 10 *® 

Pressure, terr 

f'ig. 2. Thermal dlffusivity of sample 77017,24 as a function of ombfeni gas (air) pressure, 
The data arc taken at lemperatures around 400*K. 


pressure is shown in Fig. 2. As e,xperimentally .jhown by Fujii and Osako (1973) 
for lunar rocks and also shown by Wechsicr and Glaser (1965) for particulate 
terrestrial rocks, the thermal diffusivity of sample 77017,24 decreases with 
pressure ranging 1 atm to 10 “* torr and stays constant at pressures below lO"’ torr. 
Similar behavior of the thermal diffusivity has been also found for sample 
70215,30 by us artti for Apollo 12 rock 12002,85 by Horai and Winkler (1974). 
These observatiotis mdicate that a vacuum of 10”* to 10"* torr is low enough to 
simulate the lunar environment, as far as the thermal diffusivity is concerned. In 
Fig. 3 is Shown the temperature dependence of the thermal diffusivities of samples 
77017,24 and 70215,30 in vacuums lower than 10 "* torr. The thermal diffusivities of 
both the samples decrease with temperature up to 400*K and stay constant or 
increase slightly at temperatures from 400®K to 600"K. The temperature depen- 
dence of the thermal diffusivities can be well expressed by an equation of the 
following type, 

•fc = A+|+CT* ( 1 ) 

The constants A, B, and C arc obtained by least-square fit: 

As* 0.5235 X 10 ”’ cmVsec 
B 0.1 561 cm*' *K/scc 7^,7 ,4 

C - 0.1367 X 10”'* cm’/sec-‘K* 
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CALCULATIONS Of THERMAL DlfrUSIVITY AS A 


ruNCTXOM or temperature in lunar 
•ASALT IN RANOE lOO^^K * 2000^K ♦ 


M? 97 Proari 


Th«nMl 

_ Diffusivity 

ISa C.2/..C 


k • A B . ♦ CT 




P * • 

• • • • 


1 1 «aV. 




. • . . *v • ♦ • ^ 

• • • • 


. • • • 


2»ooV. 


A ■ O.S23S X 10 ^ cm^/*«c 


B • 0.1S61 


2 

cm K/s«c 


C ■ 0.1367 X 10 cm^/s«c 


Ragist«r 0 
R«glst«r 1 


Register 2 


M«thod» Enter “■enperatiur* and key A and program prints the diffuslvity using 


aquation (1) 





A-30 


APPENDIX 3.3 

HAND CALCULATOR PROQRAM 
TO DETERMINE TEMPERATURE FIELD 
IN A LAMINAR PLOW OF 

constant viscosity with 

UNIFORM HEAT INPUT AT THE 
BASE OF THE FLOW 


o 
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APPENDIX 3.4 


Co«PTi«on of condition! for I—dnT and turbul«nt flow in f vral 
tluid« ot q«ologic«l intT«»t« 


B4s1c Physical PT*ii>»t«r8 


Vl«cosity co«fflcl»nt 


air • 0.14 eg* units 

■ 140 9 /(cm*(sac) 

poiss ■ cgs unit of absoluts viscosity ■ g /sac x 

Poisa ■ cgs unit of abs. vise. ■ gm/sac x cm 

■ cgs unit of kin. vise. ■ g /sac x cm x dans i P) 




Viscosity of watar at 


20^C ■ 0.01 ,’’/oisa 



Flow of liquid through a tuba 


If tha tangantial forca exartad by a layar of fluid upon ona adjacant layar 
is ona dyna for a spacti rata of variation of tangantial valocity than viscosity 
is ona poisa 


2 

cm 


• -^) ■ dyna-saconds • gm cm 

sac cm j * 1 — 

cm sac 


-6 o 

Air “ 182 X 10 poisa at 18 C 

Watar “ 1.002 x 10 ^ poisa at 2C^C 
Lava ■ 6. S to 7.5 X 10^ poisa at 1130 to 113S^C 


sac “ gm 

2 

cm cm sac 


P F.56 Handbood of Chandstiy ind Physics 
P F.49 Handbook of C3iamistry and Physics 

Moora, H.J. and Schabar. G.C. (197j) , An Estimata of tha Yield Strength of 
the Imbriun Flows. Proc. Lun. Sci. Conf. 6th » p. 101-118. (P. 105 quoting 

work by Shreve et al . , 1968). 



APPENDIX 3.4 

Compariaon of condibioni for laminar and turbulent flow in s#vral 
fluids of gaoloqical Intarasta 


Basic Physical Parawatars 


Viscosity coofficiant 

air ■0.14 cgs v^nits 

■ 140 g/(cm| (sec) 


jpoise ■ cgs unit of absolute viscosity ■ g /sec x 

Poise ■ cgs unit of abs , vise . ■ gm/sec x cm 

■ cgs unit of kin. vise. ■ g /sec x cm x dens (®P) 


Viscosity of water at 20”c ■ 0.01 poise 


TTprl 

8V^ 


Flow of liquid through .a tube 


If the tangential force exerted by a layer of fluid upon one adjacent layer 
is one dyne for a space rate of variation of tangential velocity then viscosity 
is one poise 


* • 

2 

cm 

cm 1 

(— ) ■ dyne-seconds ■ 

sec cm ^ 

cm 

gm cm 
2 

sec 

Air 

» 182 x 10 ^ poise at is'^c 


Water 

■ 1.002 X IP ^ poise at 20^0 


Lava 

■ 6.5 to 7.5 X 10^ poise at 1130 to 

1135®C 

^ P F.56 

^ ^ ^ fi 

Handbood of CSiemistiy nnd Physics 

- 

^ P F.49 

Handbook of Chemistry and Physics 



sec 


cm 




cm sec 


Moore* H.J. and Schaber, G.C. (1973) * An Estimate of the Yield Strength of 
the Iinbrium Flows. Proc. Lun. Sci. Conf. 6th , p. 101-118. (P. 105 quoting 

work by Shreve et al., 1968) . 




pJ3i- 

R (Rsynolds nuinbac ) » ^ ^ 

Flow i« laminar genarally fpr R^ 1000 
Plow ia turbulont genarally for R 1000 


Substance 

/ 3 

o/cm 

Poise. 


L 

cn 

V . 

cm/sec 

R 

Air (earth) 

'l 

3^ 

1.205 X 10^ 

1.82 X 10“® 

6.59 

* 

10^ 

10 

10 

10 

6590 

659 

Water 

1.0 

, 1.002xlU^ 

IQO 


10 

1 

1 

1000// 





100 

10 

lofi 

Lava 

' ' 1 

3.0 

1 X 10^ 

3 X 10 

-3 

10^ 

3 X 10^ 

' 

900 

Air (Mars) 

1 X 10“® 

160 X lO"® 

6.2 X 

10*^ 

10^ 

3 X 10^ 

1 

1.8 X 


Handbook of Chem. Rhys. P F-ll 

Carr, M.H. , The Role of Lava Erosion in the Foinnation of Lunar 
Rilles and Martian Chauinels, Icarus, 1-23, 1973. 

Viscosity is indep of pressure but does depend on tenit CaT ) . Assume 
230°K in atmosphere on Mars. Yovorsky and Detlaf (sea Handbook of 
Physics, Mir Publishars, Moscow). 


o ' " ' 
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Th* conclusion from this comparison is that water and air flow will ba 
turbulent in all conditions of the scale of those likely to be ;nd in nature 
whereas there is a regime of relatively small scale or low velocity flow 
or both where lava will flow in a laminar fashion. Because of the lower 
density of the atmosphere on Mars* it is much closar to being in a laminar 
regime than on the earth • the atmsopheric viscosity is about the same. 
However* it still isn't in the same ball park. Interestingly* the kinematic 
viscosity is approaching that of lava. 





APPENDIX 3.5 


COMPUTER PROGRAM FOR COMPUTING 
EQUILIBRIUM SLOPE PROFILE FOR SIMPLE 
PHENOMENOLOGICAL MODELS OF 
LAVA EROSION 


This is a program for a Hewlett Packard 97 programmable 
calculator that determines the slope profile for a slope that 
has reached an equilibrium shape under the action of thermo- 
mechanical erosion. 



* This program segment can be used In conjunction with the HP-97 Numerical 
Integration (Simpson's Rule) program to evaluate the given integral in 
terms in Po^ He, m and H. A future program should be structured to 
first calculate Po and Ho from equation (13) and (19) and store them as 
constants for use in the integration procedure rather than repetitively 
recalculating them. 



APPENDICES to TASK 4 


Abf tract of Hodgion (1969) worK on la^m flow simllltuda 
using carbowaxas with, data on physical propartias of 
carbowax and lava flows and modal paramatars for lava 


Computar simulation tast programs 




APPENDIX 4 . 1 


ABSTRACT OF HODGSON (1969) WORK ON LAVA FLOW SIMILITUDE USING 
CARBOWAXES WITH 3kTk ON PHYSICAL PROPERTIES OF CARBOWAX AMD 
LAVA FLOWS AND MODEL PARAMETER VALUES FOR WAX AND LAVA 


ABSTRACT I 

Th9 mechanism of emplacement of viscous flows was invest!- | 

gated experimentally with Carbowax materials. The design | 

requirements for modeling the empiacement of nature! lau^a flows | 

of high and low SiOjj content were considered. These requirements 
were found to be very complex and the laboratory wax models 
contained distortions. However, Carbowax flows are affected by 
similar heat and mass transfer processes as natural flows and 
behave in similar ways. Hence, they provide a qualitative and 
possibly semi-quantitative insight into the development of 
natural lava flows. 

The models of emplacement and many structural features 
observed in the experimental flows are similar to those in 
nature. The influeence of volume, slope, initial velocity, and 
extrusion rate on the Carbowax flows was investigated. The 
following scaling relationship was developed for maximum flow 
length relative to any reference flow as a function of the physical 
parameters investigated. j 



m ' ■ 

B * 0.382 * .136 • Q « 

b;» 0.4T9 > .085 (j a 0.459 ± .014, 

. * -0.174 * .039 


The positive value of d would apply to extrusion rates in nature | 
of less than 430 million m^/day. The negative value applies to 
higher extru,sion rates . 
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APPENDIX ^ (cont»d.) 

(Ftom wlddicioinb* «nd McGatehin* 2.97E) 


Model ParaBjeter Values for Lavi 
and Carbowax Flows 


PI ’ 
Terms 


Basic 
Lava Flows 

Acidic 
Lava Flows 

Carbovvax 
4000 Flows 

Carbowax 
20M Flows 

Wfdth 

Hatio 

b 

t 

.20 

,20 

.20 

* 

.20 

Cepth 

Ratio 

d 

X 

.02 

.10 

.00 

• 

.20 

Roughness 

r 

variable 

variable 

sinooth 

smooth • 

Froude 

number 

v2 

gii 

10-3 

io-« • 

10"^ 

10’® 

Reynolds 

Number 

gvd 

V 

10-’ 

10"* • 

10-^ 

• 10-* 

Weber 

Number 

£v!i 

c 

10* 

lo’ 

10“^ 

10-® 

Nusselt 

Number 

•>0* . 
ir 

l(10*)hj 

l(10*)he 

6(103)hg 

®0l)3)l>e ■ 

Prandtl 

Numbor 

^ • 

. ,19* 

loll 

10^ 

10* _ 

Eckert 

Number 

rj" 

p 

lO'® 

10-19 

10“® 

10-11 
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APPENDIX 4.2 

FINITE ELEMENT MODELLING EXERCISES 




During the time period from Aug. to Nov. 19 79# several pro- 
grams were written and experimented with to develop an under- 
standing of two dimensional finite element models. They fall into 
two categories: (pose related to the vorticity transport 

nT .i. «. 2 — 

equation — ^ » v A z, and those which solve the Poisson 
equation « 7* The overbars represent dimensional quantities. 

A complete solution to the flow problem involves successively 

// • 

Solving each of the two equations for each time stepe 

A basic finite difference form of the vorticity transport 
equation is (see Roche, Equation 3.165): 

^ .n _n _n ^ _n ^_n 

^ H ^ .. (‘^ 


n+1 

*'i - ^i 


.. (*^i+i - 4-1) . ..(4+1 * 4-1 - ^'^1) 

-u y; ’ .: + a — i. — .i — 

(Ax)^ 


Tax 

ADVECTION TERM 


4.2.1 


DIFFUSION TERM 


where u * velocity 

ct » diffusion constant 
superscripts refer to the time step 
subscripts refer to the cell number 
this is a centered time, centered space (Ctcs) equation, since 
the average in the 'time' term (LHS) is centered at the present 
step, and the average in the two 'space* terms (RHS) is centered 
at the cell whose new value is being computed. 

This method is unconditionally unstable (Roche, p. 36) . 
However, .if only the advection term is considered, the method 
yields stable solutions. To get back the diffusion term, 
the method of Du Port and Fankl was used. In this method, 
the diffusion term is replaced with the expression 
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II 


0 . 1 


n-hl ^ ^i-1 1^1 * ^L ) 

(Axp 


4.2.2 


Th# only dlfferdnce is in tho lAst; term in the numerator, which 
becomes the sum of past and future values at the cell in question. 

This is the Du-E*ort Framhl leapfrog method, and the pro- 
grams employing it are DFFLF (the one dimensional case) and 
2DFFLF (the two dimensional case) « The programs as written 
contain the boundary and initial conditions, and also contain 
the condition u » 0. Thus they, at present, deal only with 
diffusion. In fact, the 2D progrzun does not contain the 
advection term at all. The programs were tested (and compared 
with analytic solutions graphed Oy the programs ANATST and 
BNATST) with two sets of initial Conditions: c » 100, and 

II 

The Du-Fort Franlcl leapfrog method is not transportive 
and introduces phase erros^s when u * 0, so in the next 
exercise, the Du-Fort Frank! programs were modified to in- 
corporate second upwind differencing in the advection terms. 

This is essentially a forward space method, where the direction 
is decided by the sign of u. The equation (for advection only) 


IS : 


.n+1 


- ? 


n 




4.2.3 


At Ax 

Where v^ is the average of the velocity at cell i and cell i t 1 





V. is the average of the velocity at i and cell i •f 1 
is the average of c;” and 

is the average of c” and 

Note that this method is a forward time method. This 
maizes it incompatible with the Du-Fort Frankel. Nonetheless, 
hybrid programs were written by replacing the advection term ^ 
of equation 4.2.1 with the RHS of equation 4.2.3. These programs 
are called SUDLP and 2SUDLF . One of the problem^ii with these 
programs is that the advection terra updates itself only every 
other iteration, due to the leapfrogging inherent in Du-Port 
Frankel, and effectively proceeds at half speed, while the 
diffusion term proceeds at full speed. 


The other progreuns solve the Poisson equation 
emd use the method of successive over* relaxation. This is 
an iterative method which converges asymraetrically on the 
solution. The finite difference equation is: 


fk+1 

ij 


+ 

ij ' 


fa) 

2 ( 1 + 6 ^ 


tU.j * 




3-3 


- 2 


4.2.4 
2 ) 




0 ~ » cell size aspect ratio 
0 ) s relaxation parameter 1 ^ to < 2 


where 
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U 


(t) controls the speed of convergence, and the optimum u to use 
depends on the geometry and boundary conditions of the problem, 
and can be found analytically only in a limited number of eases , 
For the rest of the cases, the opt:lm'3m u> must be found by trial 
and error. 

Two programs were rewritten that solve Poisson's equation 
with Dirichlet conditions. They are SOR and S0RNQ3. The di£> 
ference between them is the boundary conditions. SOR has a 
source in the center, and has the outside held at zero potential, 
while S0RNQ3 has no sources, and has two facing boundaries 
at a potential of 100, and the other two boundaries held at 
zero. The output of these two programs was compaired to the 
analytic solution as plotted by SORTST and S0R2T respectively. 

A third version of the program, NEUSOR, solves Poissons 
equation with mixed boundasry conditions and no sources. The 
boundary conditions are Neuman on two facing edges, and Dirichlet 
on the other two facing edges. 

A program for finding the optimum u for the NEVSOR program 
was also written. It employs a decimal search, and counts the 
number of iterations required for convergence of a particular 
problem at each u. A binary search was not employed because 
if it were, the program would spend a large part of it's time 
in overrun conditions. 

Experience gained in exercising these programs was used 
in re- formulating the lava erosion problem, and writing the 
programs needed to solve it. 


LIST OF PROGRAMS 


DFPLP 

SUDLF 

AMATST 

BNATST 

2DFFLF 

2SUDLF 

SOR 

SORNQ3 

SORTST 

SOR2T 

NEUSOR 


solves 1-D vortlcity transport equation using Du-FOrt 
FranJcel leapfrog method. 

solves 1-D vorticity transport equation using Second 
Upwind Differencing on the advection term, and the 
Du-FOrt. 

graphs analytic solution of heat flow equation with 

^ ■ *1“ issltir 

graphs analytic solution of heat flow equation with 
IC's T ■ 100, except T * 0 at boundary s. 


solves 2-D vorticity equation using Du-Fort Frankel 
leapfrog method, but with no advection term in 
equation. 

solves 2-D vorticity transport equation using second 
upwind differencing on the advection term, and the 
Du-Fort Frankel method on the diffusion tenn. 


solves Poisson's equation with boundaries held at 
zero potential and a source in the center, using the 
Successive Over-relaxation. 


c ■ . 

graphs analytic solution of the problem solved by SOR. 

graphs analytic solution of thva problem solved by 
SORNQ3. 

solves Posson's equation using .Successive Over- 
relaxation for mixed Dlrchlet at^d Neumann B.C.'s, 
and no sources. 

Finds optimvun as to use for NEUSOR. 






NWSOR 
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0Q02 

0005 

00 o4 

0005 

•0006 

0007 

0000 

0000 

0010 

0011 

0012 

OOlS 

ooia 

0015 

0016 
0017 
OOlS 
0010 
0020 
0021 
0022 
0023 
0020 

0025 

0026 
0027 
0026 
0020 

0030 

0031 

0032 

0033 
0030 

0035 

0036 


0037 
0036 
0030 
004 0 

0041 

0042 
004J 

0044 

0045 

0046 


APPiPUDIX 4.3, FORTRAN Program for two dimansional modal of 
thermal aroslon for a plane laminar constant viscosity flo^l with 
a yield temperature. 

functions inline 4REl ALL 

INTECER pP»00iPRINT»PT, FLOOR*0*AA4f FLIP* HC8 

REAL MiXfHELTEMfLENGtH, KAPPA 
e ALL OTHER VARIA8LE8 FOLLOW XJKLHN CONVCNTION 

DIMENSION T(12f l7f2)fPT(12«17)90PT(i7*l2lfF(12«17)fV(i7)fDV(11 
lfV(17) 

10«5 

AAA«6 

THRB1200. 

Mix«o.u*oa 

AN6S30* 

THs(ANCy 180. 163,141592654 

P»s30 

TMAX»96, 

ROCaO,25 

VMAXbIOQO. o 

OT8400. 

JM«17 

JMMUJMwI 

ilHH2«JH»2 

J««3aJM«3 

AL9l.3S6E»12 

PRINT** I 

ITERB^l 

tIMC«-l,6DT 

MELTE^^biSOO. 

A^^8TEMB^0, 

FLOOR* U 

IFLPlaFLOORAl 

6«980. 

M«150. 

length«o.ie 0» 

0«1 
NS 2 

0*«LENCTH/10. 

KAPPA«0.04 

WRITE(AAA,250) kappa, ROC,G,ANG,H*VHAXiMIX 
250 pOHHAT(//40X»6HKAPpAa #Fl3.5/40X#6MRnC* »Fl3.5/ 

J 40Xf6H6* ,F13,5/40X.6MANC« ♦F13,5/40X> 

2 6HH ♦Fl3.5/40X,6HyMAX* iFlS.S/ttOX* 

3 6HMIX* ♦F13,5//1 

00 1 !*1 * 12 

00 2 J«1 fFLOnR 
T(I*JfO)*AMBTEM 
2 CONTINUE 

00 3 JsIFLPl ♦ JM 
T ( I ♦ J»0)*MELtEM 

? CoUtImIIe ORIGMAt PAQE I? 

REAO(2,101) fOV(J) , J«2, JMMJ 1 OF POOR Qy/\Ui » 

101 FORMAT C15F«.n 




I8N nonf 00 1050 


ISN 

nous 

1050 

Dy(J)«t50,*DV(J) 


ISN 

00/*9 


Y<jM-n,uvUM-i)/a, 


ISN 

0050 


DC 5} JJ*2*JHM2 


ISN 

005} 




ISN 

0052 

51 

y(j)«v(jAn*ov(j*n/2,*oYtJ)/2* 


ISN 

005J 


on 56 Js2vJMH1 


15N 

0050 


17 (y(J)-H) 57.58.58 


ISN 

0055 

58 

V<J)>0 


IIN 

0056 


GOTO 56 


UN 

0057. 

57 

VIJ)sVHAK6(},ii(Y(J)/H)9*2) 

1 

IIN 

0050 

56 

continue 

ft 

! 

ISN 

0059 

430 

iTENBlTERY} 

ISN 

0060 


timc«txhe*ot 

1 

ISN 406} 


DO 5 JB2trL00R 

1 

1 

ISN 

0062 


T(}»J*0)iT(2*J»0) 


ISN 

0063 
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PROPOSAL SUMMARY 


PRINCIPAL INVESTIGATOR; 
(Name » address, tel«no.) 


James A» Cutts 

Planetary science inst/scxence Applicatioi 

^ 283 S. Lake Ave. Suite M 

ii Pasadena, cA silAl 


CO-INVKSTIGATORS: Karl R^ Blasius 

(Name only) Ciar)c R, cnaproan 

Wm. James Roberts 

Title: Geophysical Constraints on Lunar t Planetary Volcanism 
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ABSTRACT (Single- Spaced, type within box below. Include: 
a. Brief statement of the overall objectives and justification 
of the work; b. Brief statement of the accomplishments of the 
prior year, or "new proposal"; c. Brief listing of what will be 
done this year as well as how and why; and d. One or two of your 
recent publications relevant to the proposed work*) : 


al We propose a one-year follow on to our investigations of 
the role of volcanism in forming features of planetary sur- ’ 
faces . Its goals are to refine and extend the qualitative 
and quantitative understanding of the formation of volcanic 
features. Task 1 is a continuation of the study of the 
physical mechanism of thermal lava erosion with, application 
to the origins of some lunar sinuous’ rilles and martian 
channels. Task 2 is a continuation of an investigation of 
the origin of small lunar' craters and its implications for 
mare basalt petrogenesis . 

b) . Through March 1979 we performed a classification of 
small crater features, catalogued central volcanic constructs, 
and gathered basic physical data on martian central volcanic 
constructs. 

c) In Task 1, physical modeling studies will be used to test 
the validity of thermal erosion codes. A solid-liquid wax 
system will be used. Initial phenomenological investigations 
of flow responses to obstacles, channel sinuosity, and 
variations in substrate characteristics will also be conduct- 
ed. In Task 2 data gathered in previous years will be pre- 
pared for publication and a simple numerical model of small 
crater production and destruction will be used in data 
interpretation . 

d) Clark R. Chapman, Jayne C. AUbele, Wra. James Roberts 
and James A. Cutts. Sub-kilometer lunar craters; Origins, 
ages, processes of degradation and implications for mare 
basalt petrogenesis. Lunar and Planetary Sciences Conference 
X Extended abstracts, 1979. 



